05 - Raster manipulation
Instructions
- Activate
rasterandrgdallibrary - Use
GDALinfo()to inspect the properties of the Aster image raster in the data folder (the file name is “Aster.tif”). What can you learn from this inspection about its bands, resolution, and projection? - Read in the Aster image. Review Week 02 raster loading if unsure about how.
- Use the
inMemory()function on the elevation object to determine if the data has been read in. - Use the
head()function to look at the first few values from the elevation raster. - Use the
getValues()function on the elevation object to read in all the data. - Use the
hist()function to create a quick histogram of the elevation values. Note the pile of values near -9999, these should beNA(any idea why?) and we will address this later.
# Library
library(raster)
library(rgdal)
# Read in the elevation layer
elevation <- raster("C:/Users/alexp/Documents/School/Spatial_analytics/HW/HW3/cds-spatial/data/Aster.tif")
# Check if the data is in memory
inMemory(elevation)[1] FALSE
# Use head() to peak at the first few records
head(elevation) 1 2 3 4 5 6 7 8 9 10 11 12
1 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
2 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
3 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
13 14 15 16 17 18 19 20
1 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
2 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
3 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
[ reached getOption("max.print") -- omitted 7 rows ]
# Use getValues() to read the values into a vector
vals <- getValues(elevation)
# Use hist() to create a histogram of the values
hist(vals)# all I could find is that 'the Values of -9999 are assigned to missing or
# cloudy data' taken from
# https://lpdaac.usgs.gov/documents/120/ASTERGED_User_Guide_V3.pdf also '-9999
# for void pixels and 0 for sea-level water body' from
# https://lpdaac.usgs.gov/documents/434/ASTGTM_User_Guide_V3.pdfInstructions
- Check that the package
rasterand the objectelevationare loaded in your workspace. - Plot the
elevationraster usingplot(). - Set up a three-column matrix with the
cbind()function and values -10000, 0, NA. - Use the matrix and
reclassify()to reclassify values below 0 toNA. You will need to use the argumentrcl. - Plot the reclassified elevation layer to confirm there are no values below 0.
# Plot the elevation layer to see the legend values
plot(elevation)# Set up the matrix
vals <- cbind(-10000, 0, NA)
# Reclassify
elevation_reclass <- reclassify(elevation, rcl = vals)
# Plot again and confirm that the legend range is 0 - 2400
plot(elevation_reclass)Instructions I
- Load
sflibrary - Create
moundsobject from the shapefile “KAZ_mounds.shp”. For a refresher on vector loading, check out Week 02 instructions. - Verify that
moundsCRS matches theelevationand fix withst_transform()if not. - Create a bounding box around the mounds with
st_make_grid()following Week 03 guidelines. - Plot the
mounds_bband themoundsto see how these two objects relate. - Crop the
elevationlayer by themoundsobject and create a smallerelevobject - Plot the
elevandmoundsand themounds_bbtogether.
# Load `sf` library
library(sf)
# Create `mounds` object from the shapefile
mounds <- st_read("C:/Users/alexp/Documents/School/Spatial_analytics/HW/HW3/cds-spatial/data/KAZ_mounds.shp")Reading layer `KAZ_mounds' from data source
`C:\Users\alexp\Documents\School\Spatial_analytics\HW\HW3\cds-spatial\data\KAZ_mounds.shp'
using driver `ESRI Shapefile'
Simple feature collection with 773 features and 5 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 352481.3 ymin: 4712325 xmax: 371282.4 ymax: 4730029
Projected CRS: WGS 84 / UTM zone 35N
# Verify that `mounds` CRS matches the `elevation`
st_crs(mounds)Coordinate Reference System:
User input: WGS 84 / UTM zone 35N
wkt:
PROJCRS["WGS 84 / UTM zone 35N",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 35N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",27,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
ID["EPSG",32635]]
st_crs(elevation) #they matchCoordinate Reference System:
User input: WGS 84 / UTM zone 35N
wkt:
PROJCRS["WGS 84 / UTM zone 35N",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 35N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",27,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["Between 24°E and 30°E, northern hemisphere between equator and 84°N, onshore and offshore. Belarus. Bulgaria. Central African Republic. Democratic Republic of the Congo (Zaire). Egypt. Estonia. Finland. Greece. Latvia. Lesotho. Libya. Lithuania. Moldova. Norway. Poland. Romania. Russian Federation. Sudan. Svalbard. Turkey. Uganda. Ukraine."],
BBOX[0,24,84,30]],
ID["EPSG",32635]]
# Create a bounding box around the mounds
mounds_bb <- st_make_grid(mounds, n = 1)
# Plot the `mounds_bb` and the `mounds`
plot(mounds_bb)
plot(mounds, add = TRUE)# Crop the `elevation` layer by the `mounds` object and create a smaller `elev`
# object
elev <- crop(elevation, mounds)
# Plot the `elev` and `mounds` and the `mounds_bb` together
plot(elev)
plot(mounds_bb, add = TRUE)
plot(mounds, add = TRUE)Instructions II
- Create
surveyobject from shapefile called “KAZ_units.shp” - Project the
surveyobject to match the newelevraster withst_transform()if needed. - Compute the area of the survey with
st_area()and save this object asareas. What units are these? - Filter the survey units to only those above 30000 square meters with the
filter()function. You will need to wrapareasinunclass(). Save assurvey_big. Remember to have the tidyverse or dplyr library attached forfilter()to work properly. Also sf might interfere so specify the dplyr::filter() if needed.
library(dplyr)
# Read in the survey object
survey <- st_read("C:/Users/alexp/Documents/School/Spatial_analytics/HW/HW3/cds-spatial/data/KAZ_units.shp")Reading layer `KAZ_units' from data source
`C:\Users\alexp\Documents\School\Spatial_analytics\HW\HW3\cds-spatial\data\KAZ_units.shp'
using driver `ESRI Shapefile'
Simple feature collection with 7708 features and 25 fields (with 1 geometry empty)
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 352181.7 ymin: 4712923 xmax: 370892.1 ymax: 4729960
Projected CRS: WGS 84 / UTM zone 35N
# projection is not needed
# Compute the area of the survey
areas <- st_area(survey)
# assessing from the Environment menu areas are in sq meters
# Filter to survey with areas > 30000
survey_big <- dplyr::filter(survey, unclass(areas) > 30000)Instructions III
- Review the plot of
elevraster. - Plot the geometry of the
survey_bigover it. - Mask the
elevlayer withsurvey_bigand save aselevation_mask. This may take a couple of seconds. - Review the plot of
elevation_mask.
# Plot the elevation raster
plot(elev)# Plot the geometry of survey_big
plot(st_geometry(survey_big))# Mask the elev layer with survey_big and save as elevation_mask
elevation_mask <- mask(elev, mask = survey_big)
# Plot elevation_mask -- this is a raster!
plot(elevation_mask)Question:
1. What extent does the elevation raster default to after cropping by mounds vs by masking by the large polygons of survey_big? -after cropping it creates a box of an area where all the mounds are situated but after masking, while the area is the same we can only see the elevation for the specific sites
2. Why do we not use the mounds bounding box to crop the elevation raster? -I think it is not precise enough
Instructions
- Load the East half of the two IKONOS images of the Kazanlak Valley from (“KazE.tif”).
- Plot the file you read in with the appropriate function for a multi-band raster.
- Determine the raster resolution using
res()and number of raster cells in the layer withncell(). - Aggregate the IKONOS image using the default for fun and a factor of 10 and save the new raster to
east_small. - Plot the new raster with
plot()for comparison to the old version. - Determine the new raster resolution and the number of raster cells.
# Load the IKONOS raster
east <- brick("C:/Users/alexp/Documents/School/Spatial_analytics/HW/HW3/cds-spatial/data/KazE.tif")
# Plot the IKONOS raster
plot(east)# Determine the raster resolution
res(east)[1] 10 10
# Determine the number of cells
ncell(east)[1] 2374596
# Aggregate the raster
east_small <- aggregate(east, fact = 10)
# Plot the new east layer
plot(east_small)# Determine the new raster resolution
res(east_small) #100 100 vs previous 10 10[1] 100 100
# Determine the number of cells in the new raster
ncell(east_small) #23940 vs previous 2374596[1] 23940
Instructions
- Ensure mounds are still in memory.
- Use the
rasterfunctionextract()to determine the elevation at each of the mounds. Assign the extracted value into aelevcolumn in the mounds object. Beware that the extract() function exists across multiple packages, so it’s wise to use raster:: in front of it. - Look at the
moundselevation column through the plot() function as well as histogram. Do theextract()results make sense? Theelevationlayer values represent meters above sea level.
# Extract the elevation values at the sites
mounds$elev <- raster::extract(elevation, mounds)
# Look at the mounds and extraction results
plot(mounds["elev"])hist(mounds$elev) #I don't see anything wrong with the results so farInstructions
- Make sure elevation object exists in memory (“Aster.tif”, it is a single-band raster).
- Read in the prominence raster layer from “prominence.tif”; it is also a single-band raster.
- Plot the prominence object. Do the legend units make sense?
- Specify function f, where you look for elevation values over 400 msl and below 650m (mounds don’t appear above that elevation) and prominence values over 60%
- Call
overlay()onelevationandprominence. Set thefunargument tof. - If you have not cropped the two rasters to the same extent yet, you can do so inside the
overlay()function
# Check in on the elevation and read in prominence layer
elevationclass : RasterLayer
dimensions : 2484, 2592, 6438528 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 313183.7, 390943.7, 4676606, 4751126 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=35 +datum=WGS84 +units=m +no_defs
source : Aster.tif
names : Aster
values : -9999, 2347 (min, max)
prominence <- raster("C:/Users/alexp/Documents/School/Spatial_analytics/HW/HW3/cds-spatial/data/prominence1500m.tif")
# Plot prominence
plot(prominence) # supposing the units are percentages as mentioned above, yes they make sense# Function f with 2 arguments and the raster math code
f <- function(rast1, rast2) {
rast1 < 650 & rast1 > 400 & rast2 > 60
}
# Align the extent of the two rasters with crop()
elev2 <- crop(elevation, prominence)
# Do the overlay using f as fun.
elevation_prom_overlay <- overlay(elev2, prominence, fun = f)
# Plot the result (low elevation and high prominence areas)
plot(elevation_prom_overlay)elevation_prom_overlayclass : RasterLayer
dimensions : 1201, 1201, 1442401 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 340153.7, 376183.7, 4700126, 4736156 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=35 +datum=WGS84 +units=m +no_defs
source : memory
names : layer
values : 0, 1 (min, max)
Questions:
3. What is the actual value range in the prom_el_overlay raster? -the range of values is: 0, 1 (min, max)
4. What area is covered by each cell? -I don’t know if I assume correctly but if the resolution is 30 30 than the cell covers 30m x 30m as the units in coordinate system are meters
Task 7 - HOMEWORK: Inspect the most prominent mounds
In the above exercise, we produced an elevation-prominence overlay. Mounds and other sites that sit in this overlay enjoy a strategic position vis-a-vis the rest of the valley. Which ones are they, however, and what are their real prominence values? It is hard to see at the scale of the valley and it would be good to pull the sites out.
Find the mounds that enjoy the most prominent locations as well as those that feature in the elevation-prominence overlay raster. Produce a list of the ID numbers (TRAP_Code) of all the overlay mounds and 25 most prominent mounds and plot them (expressing their prominence somehow) .
Instructions
Find or recreate the
moundssf object,prominenceandprom-el-overlayrasters.Plot the
prom-el-overlayand themoundson top of each other to check visual overlap. Do the same withprominenceraster andmounds.Extract values from the elevation-prominence overlay raster and from the prominence raster at mounds and write them to two columns:
– call the first column
prom_el_overlay– name the second columnprominenceMake an object of mounds that sit within these strategic high-visibility locations. How many are there?
Make an object of 25 mounds with the highest prominence values. Which
TRAP_Codeids are included?Plot both these sets of mounds in
mapview()and compare their locations.
# Plot the `prom-el-overlay` and the `mounds` on top of each other
plot(elevation_prom_overlay)
plot(mounds, add = TRUE)# same with `prominence` raster and `mounds`
plot(prominence)
plot(mounds, add = TRUE) #I can't really see the mounds anymore, it's too cluttered# Extract values from the elevation-prominence overlay and from the prominence
# at mounds and write them to two columns
mounds$prom_el_overlay <- raster::extract(elevation_prom_overlay, mounds)
mounds$prominence <- raster::extract(prominence, mounds)
# Make an object of mounds that sit within these strategic high-visibility
# locations
best_mounds <- dplyr::filter(mounds, prom_el_overlay == 1)
# taking all mounds where the value for prom_el_overlay was 1, meaning they had
# higher than 60% prominence and desired elevation as we filtered with f before
nrow(best_mounds)[1] 132
# there is 132 mounds withing the most visible areas
# Make an object of 25 mounds with the highest prominence values
highest_prom_mounds <- mounds %>%
arrange(desc(prominence)) %>%
slice(1:25)
# writing out all the unique Trap codes included
unique(highest_prom_mounds$TRAP_Code) [1] 2002 3112 3105 4057 4058 3722 4059 3417 4056 3276 3024 1020 4087 1019 3026
[16] 2028 4066 5002 4086 1071 3275 2003 2004 3703 3420
library(mapview)
mapview(best_mounds, col.regions = "green") + mapview(highest_prom_mounds, col.regions = "red")# we can see many of these points are actually overlapping so lets make it
# easier to see the overlapping ones
# create object for overlapping points
overlap <- st_join(best_mounds, highest_prom_mounds, left = FALSE)
# add the object to mapview
mapview(best_mounds, col.regions = "green") + mapview(highest_prom_mounds, col.regions = "red") +
mapview(overlap, col.regions = "orange")Question:
5. How do the mounds with high prom_el_overlay values differ from those with high prominence? -some of the mounds with highest prominence were not high but instead located close to the lake at lower elevation
Task 8 - OPTIONAL HOMEWORK: Sum raster values across polygons - assessing prominence of past settlements
Previously you have extracted values from computed rasters at points. More often, you will need to extract values using more complex vector shapes where you will need to decide how to summarize the raster values that fall in. In this exercise, you will load the layer of ancient settlements for the Kazanlak Valley and sum the overlapping elevation_prom_overlay values within the settlement polygons. Then you will find out which of the sites contained at least 0.5 ha of strategic high-visibility area.
Instructions
- Create a
sitessf object by reading in polygons from “KAZ_nuclei.shp” - Plot the elevation-prominence overlay and the sites on top of each other
- Use
maskfunction to create asites_maskraster, where raster cells under sites’ polygons retain their original values, while cells outside the polygons are turned toNA. - Use
cropfunction to cropsites_maskraster to the extent of the sites’ polygons. - Plot the result ( I recommend
mapview()) - Produce a list of sites that have at least 0.5 ha (ie. 6 grid cells) within these strategic high-visibility locations. (Hint:you can first extract a list of overlapping sites and then use
sapply()to filter out the ones > 6)
# Read in the site polygons
sites <- ______("data/KAZ_nuclei.shp")
# Plot the sites
plot(elevation_prom_overlay);plot(__________)
# Create a raster mask from the elevation_prom_overlay using the sites
sites_mask <- _____________
# Crop the raster mask to the sites
sites_strat <- _____________
# Plot and compare the results
plot(sites_mask)
plot(sites_strat)
# Produce a list of the most prominent sitesLovely work! Manipulating rasters and mastering raster-vector interactions is the bread and butter of spatial data wrangling :)